Load packages

## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2", "kableExtra", "knitr", "formatR", "randomForest", "Rtsne", "MASS", "adehabitatHR", "sp", "GGally", "spatstat", "raster", "vegan", "brms", "lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes", "cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace")

aa <- lapply(x, function(y) {
  
  # get pakage name
  pkg <- strsplit(y, "/")[[1]]
  pkg <- pkg[length(pkg)]
  
  # check if installed, if not then install 
  if (!pkg %in% installed.packages()[,"Package"])  {

      if (grepl("/", y))  remotes::install_github(y, force = TRUE) else
    install.packages(y) 
    }

  # load package
  a <- try(require(pkg, character.only = T), silent = T)

  if (!a) remove.packages(pkg)
  })

Functions and global parameters

opts_knit$set(root.dir = "..")

# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE, tidy = TRUE)

read_excel_df <- function(...) data.frame(read_excel(...))

trait_space_plot <- function(X, dimensions, group, indices, basecex = 1, title = "", colors = c("#3E4A89FF", "#35B779FF"), point.alpha = 0.7, background.indices = NULL, pch = 1, labels = c("sub-space", "total space"), legend = TRUE) {

   total_coors <- as.ppp(as.matrix(X[, dimensions]), c(range(X[, dimensions[1]]), range(X[, dimensions[2]])))
  total_space <- raster(density.ppp(total_coors))

  xlim <- range(X[, dimensions[1]]) 
  ylim <-  range(X[, dimensions[2]])
   # ylim <- if (is.null(background.indices))
   #  range(X[, dimensions[2]]) else 
   #    c(min(X[, dimensions[2]]), c(max(X[, dimensions[2]]) + (max(X[, dimensions[2]]) - min(X[, dimensions[2]])) * 0.11))
  
  plot(x = X[, dimensions[1]], y = X[, dimensions[2]], col = "white", cex = basecex, xlab = dimensions[1], ylab= dimensions[2], xlim = xlim, ylim = ylim, cex.lab = basecex, xaxs="i", yaxs="i") 

  # add background group
  if (!is.null(background.indices)){
      
  # get points  
  Y_bg <- X[background.indices, ]
    
    if (nrow(Y_bg) > 1 & point.alpha > 0)
      points(Y_bg[, dimensions[1]], Y_bg[, dimensions[2]], col = adjustcolor(colors[2], alpha.f = point.alpha), pch = pch)
    
    if (nrow(Y_bg) > 4){
  coors <- as.ppp(as.matrix(Y_bg[, dimensions]), c(range(Y_bg[, dimensions[1]]), range(Y_bg[, dimensions[2]])))
      
 raster_dens <- raster(density.ppp(coors))

  palpre <-  colorRampPalette(c("white", colors[2]))    
  colsn <- 10
  palpre2 <- sapply(1:colsn, function(x) adjustcolor(col =  palpre(colsn)[x], alpha = seq(0.1, 0.9, length.out = 20)[x]))
  
  image(raster_dens, add = TRUE, col = palpre2)
  }
  
  }
  
  Y <- X[indices, ]
  
  if (nrow(Y) > 1 & point.alpha > 0)
      points(Y[, dimensions[1]], Y[, dimensions[2]], col = adjustcolor(colors[1], alpha.f = point.alpha), pch = pch)
   
    if (nrow(Y) > 4){
  coors <- as.ppp(as.matrix(Y[, dimensions]), c(range(Y[, dimensions[1]]), range(Y[, dimensions[2]])))
      
 raster_dens <- raster(density.ppp(coors))

  palpre <-  colorRampPalette(c("white", colors[1]))    
  colsn <- 10
  palpre2 <- sapply(1: colsn, function(x) adjustcolor(col =  palpre(colsn)[x], alpha = seq(0.1, 0.9, length.out = 20)[x]))
  
  image(raster_dens, add = TRUE, col = palpre2)
    }
  
   usr <- par()$usr
  
  # legend  
  if (!is.null(background.indices) & legend)
  # legend(y = rep(usr[4] + ((usr[4] - usr[3]) * 0.05), 2), x = c(usr[1] + (usr[2] - usr[1]) * 0.25, usr[1] + (usr[2] - usr[1]) * 0.75), col = colors, legend = labels, pch = c(pch, pch), cex = basecex * 1.5, bty = "n", horiz = TRUE, xjust = 0)
  legend("topright", col = colors, legend = labels, pch = c(pch, pch), cex = basecex * 1.5, xjust = 0, bg = adjustcolor("white", alpha.f = 0.5))
  
  title(title, cex.main = basecex * 1.2)
}

# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")


cols <- viridis(10, alpha = 0.7)

 summary_brm_model <- function(x){
  
   summ <- summary(x)$fixed
  fit <- x$fit  
  betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE)  
  hdis <- t(sapply(betas, function(y)   hdi(fit@sim$samples[[1]][[y]]))
)
  summ_table <- data.frame(summ, hdis, iterations = attr(fit, "stan_args")[[1]]$iter, chains = length(attr(fit, "stan_args")))
  summ_table <- summ_table[rownames(summ_table) != "Intercept", c("Estimate", "Rhat", "Bulk_ESS", "l.95..CI", "u.95..CI", "iterations", "chains")]
  
  summ_table <- as.data.frame(summ_table)  
  summ_table$Rhat <- round(summ_table$Rhat, digits = 3)  
  summ_table$CI_low <- round(unlist(summ_table$l.95..CI), digits = 3)  
  summ_table$CI_high <- round(unlist(summ_table$u.95..CI), digits = 3)  
  summ_table$l.95..CI <- summ_table$u.95..CI <- NULL
  
    out <- lapply(betas, function(y)  data.frame(variable = y, value = sort(fit@sim$samples[[1]][[y]], decreasing = FALSE)))
  
  posteriors <- do.call(rbind, out)
  posteriors <- posteriors[posteriors$variable != "b_Intercept", ]
  
   gg <- ggplot(data = posteriors, aes(y = variable, x = value, fill = stat(quantile))) + 
  geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = HDInterval::hdi,  vline_linetype = 2) + 
  theme_classic() + 
     xlim(c(min(summ_table[ , "CI_low"]), max(summ_table[ , "CI_high"]))) +
  geom_vline(xintercept = 0, col = "red") + 
  scale_fill_manual(values = c("transparent", "lightblue", "transparent"), guide = "none") + ggtitle(x$formula)

  
  summ_table$Rhat <- ifelse(summ_table$Rhat > 1.05, cell_spec(summ_table$Rhat, "html", color ="white", background = "red", bold = TRUE,  font_size = 12),  cell_spec(summ_table$Rhat, "html"))
  
  signif <- summ_table[,"CI_low"] * summ_table[,"CI_high"] > 0
  
  df1 <- kbl(summ_table, row.names = TRUE, escape = FALSE, format = "html", digits = 3)
    
  df1 <- row_spec(kable_input = df1,row =  which(summ_table$CI_low * summ_table$CI_high > 0), background = adjustcolor(cols[9], alpha.f = 0.3))

  df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12)
  
  print(df1)
  
  print(gg)

  }
# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")

# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd", "freq.median",
    "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
    "mindom", "maxdom", "dfrange", "modindx", "meanpeakf")]

sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])

sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year, sep = "-")

sels$date <- as.Date(sels$date, format = "%d-%b-%Y")

Counts per individual

df <- as.data.frame(table(sels$ID))

names(df) <- c("ID", "Sample_size")

df <- df[order(df$Sample_size, decreasing = FALSE), ]

kb <- kable(df, row.names = FALSE)

kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

print(kb)
ID Sample_size
132YGMM 6
125YGMM 12
160YGHM 16
310YGHM 21
393YGHM 25
303YGHM 37
398WBHM 41
365YGHM 44
399YGLM 46
300YGHM 47
270YGHM 64
407YGHM 97
386WBHM 100
377WWLM 108
367WWNM 119
371YYLM 149
397YGHM 175
378YYLM 195
404WBHM 196
258YGHM 223
389WWLM 230
262YGHM 306
400YGHM 360
388YGLM 377
327YYLM 404
396YBHM 455
403WBHM 566
356WBHM 761
361YGHM 767
402YGHM 849
395WBHM 854
360YGHM 900
390WBHM 935
405YBHM 975
385YBHM 1256
394YBHM 1693

Add metadata

metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")

# head(metadat)

sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"

# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID, sels$ID)

sels$treatment <- sapply(1:nrow(sels), function(x) {

    metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]

})


sels$treatment.room <- sapply(1:nrow(sels), function(x) {

    metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]

})


sels$round <- sapply(1:nrow(sels), function(x) {

    metadat$Round[metadat$Bird.ID == sels$ID[x]][1]

})

sels$source.room <- sapply(1:nrow(sels), function(x) {

    metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]

})

sels$record.group <- sapply(1:nrow(sels), function(x) {

    metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]

})


# add week
out <- lapply(unique(sels$round), function(x) {

    Y <- sels[sels$round == x, ]

    min_date <- min(Y$date)
    week_limits <- min_date + seq(0, 100, by = 7)

    Y$week <- NA
    for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i - 1] & Y$date <
        week_limits[i]] <- i - 1

    return(Y)
})

sels <- do.call(rbind, out)


sels$cort.baseline <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


sels$cort.stress <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})

sels$stress.response <- sels$cort.stress - sels$cort.baseline

sels$weight <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


sels$breath.count <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


# remove week 5
sels <- sels[sels$week != 5, ]

Exploratory graph

Physiological parameters

breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count", "D14.Breath.Count",
    "D21.Breath.Count", "D28.Breath.Count")])

weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.", "D14.Bird.Weight..g.",
    "D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])

cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress", "D14.CORT.Stress",
    "D21.CORT.Stress", "D28.CORT.Stress")])

cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline", "D14.CORT.Baseline",
    "D21.CORT.Baseline", "D28.CORT.Baseline")])


stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round", "Treatment.Room")],
    week = breath.count$ind, breath.count = breath.count$values, weight = weight$values,
    cort.stress = cort.stress$values, cort.baseline = cort.baseline$values, stress.response = cort.stress$values -
        cort.baseline$values)

# head(stress)

stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."), "[[", 1),
    levels = c("D3", "D7", "D14", "D21", "D28"))

stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)

stress$treatment <- factor(stress$Treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

# remove 5th week
stress <- stress[stress$week != "D28", ]


stress_l <- lapply(stress$Bird.ID, function(x) {
    X <- stress[stress$Bird.ID == x, ]

    X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
    X$weight <- X$weight - X$weight[X$week == "D3"]
    X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
    X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week == "D3"]
    X$stress.response <- X$stress.response - X$stress.response[X$week == "D3"]

    return(X)
})

stress <- do.call(rbind, stress_l)

ggplot(data = stress, aes(x = day, y = weight, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Weight (g)", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = breath.count, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Breath count", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = cort.baseline, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Baseline CORT", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = cort.stress, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Stress CORT", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = stress.response, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Stress response", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

# stress_l <- lapply(unique(stress$Bird.ID), function(x){ X <-
# stress[stress$Bird.ID == x, ] X$cort.baseline.change <- X$cort.baseline - X$
# })

standard_error <- function(x) sd(x)/sqrt(length(x))

agg_stress <- aggregate(cort.baseline ~ week + treatment, stress, mean)
agg_stress$se <- aggregate(cort.baseline ~ week + treatment, stress, standard_error)[,
    3]
agg_stress$Week <- 1:4

ggplot(data = agg_stress, aes(x = Week, y = cort.baseline, fill = treatment)) + geom_bar(stat = "identity") +
    geom_errorbar(aes(ymin = cort.baseline - se, ymax = cort.baseline + se), width = 0.1) +
    # geom_line(size = 1.2) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") +
    labs(y = "Baseline CORT", x = "Week") + theme_classic(base_size = 30) + theme(legend.position = "none")

FOXP2

foxp2 <- read.csv("./data/raw/stress_IHC_counts.csv")

foxp2$Treatment[grep("HIGH", foxp2$Treatment)] <- "High stress"
foxp2$Treatment[grep("Control", foxp2$Treatment, ignore.case = T)] <- "Control"
foxp2$Treatment[grep("MED", foxp2$Treatment)] <- "Medium stress"

foxp2$Treatment <- factor(foxp2$Treatment, levels = c("Control", "Medium stress",
    "High stress"))

ggplot(data = foxp2, aes(x = Treatment, y = FoxP2.Counts.MMSt.Striatum, fill = Treatment)) +
    geom_boxplot() + geom_point() + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
    labs(y = "% FoxP2 expression", x = "Treatment") + theme_classic(base_size = 30) +
    # scale_x_discrete( labels = c(1:4)) +
theme(legend.position = "none")

Behavioral parameters

Call counts per treatment

call_count <- aggregate(selec ~ week + treatment + round, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (unique(call_count$week))[5:1])

call_count$round <- paste("Round", call_count$round)

call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

call_count$stress.response <- aggregate(selec ~ week + treatment + round, data = sels,
    mean)[, 4]

ggplot(data = call_count, aes(x = treatment, y = selec, fill = Week)) + geom_bar(stat = "identity") +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + coord_flip() + facet_wrap(~round,
    ncol = 2) + labs(y = "Number of calls", x = "Treatment") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.15))

call_count2 <- aggregate(selec ~ week + treatment + round + ID, data = sels, length)
agg_count <- aggregate(selec ~ week + treatment, call_count2, mean)
agg_count$se <- aggregate(selec ~ week + treatment, call_count2, standard_error)[,
    3]

agg_count$treatment <- factor(agg_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

ggplot(data = agg_count, aes(x = week, y = selec, fill = treatment)) + geom_bar(stat = "identity") +
    geom_errorbar(aes(ymin = selec - se, ymax = selec + se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Number of calls",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = "none")

Call counts per bird and treatment

# Call counts per treatment by individual (in color)

call_count <- aggregate(selec ~ week + treatment + round + ID, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (1:5))

call_count$round <- paste("Round", call_count$round)

call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

ggplot(data = call_count, aes(x = week, y = selec, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Call count", x = "Week") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.25))

Acoustic parameters

call_week_params <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
    time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
    modindx, meanpeakf) ~ week + treatment + round, data = sels, mean)

call_week_params_sd <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
    time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
    modindx, meanpeakf) ~ week + treatment + round, data = sels, sd)

names(call_week_params_sd) <- names(call_week_params) <- c("week", "treatment", "round",
    "Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
    "Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
    "Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
    "Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
    "Modulation_index", "Mean_peak_frequency_(kHz)")


call_week_params_sd <- call_week_params_sd[, c("Duration", "Mean_frequency_(kHz)",
    "SD_frequency_(kHz)", "Median_frequency_(kHz)", "Interquantile_frequency_range_(kHz)",
    "Interquantile_time_range_(s)", "Skewness", "Kurtosis", "Spectral_entropy", "Time_entropy",
    "Overall entropy", "Mean_dominant_frequency_(kHz)", "Minimum_dominant_frequency_(kHz)",
    "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)", "Modulation_index",
    "Mean_peak_frequency_(kHz)")]

names(call_week_params_sd) <- paste(names(call_week_params_sd), "sd", sep = ".")


call_week_params <- cbind(call_week_params, call_week_params_sd)
call_week_params$round <- paste("Round", call_week_params$round)
call_week_params$treatment <- factor(call_week_params$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))


pd <- position_dodge(0.3)

ggs <- lapply(c("Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
    "Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
    "Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
    "Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
    "Modulation_index", "Mean_peak_frequency_(kHz)"), function(i) {
    call_week_params$var <- call_week_params[, i]
    call_week_params$var.min <- call_week_params$var - call_week_params[, paste(i,
        "sd", sep = ".")]
    call_week_params$var.max <- call_week_params$var + call_week_params[, paste(i,
        "sd", sep = ".")]

    ggplot(call_week_params, aes(x = week, y = var, col = treatment)) + geom_point(size = 4,
        position = pd) + geom_errorbar(aes(ymin = var.min, ymax = var.max), width = 0.3,
        position = pd, size = 1.7) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
        facet_wrap(~round, ncol = 2) + labs(y = i, x = "Week") + theme_classic(base_size = 30) +
        theme(legend.position = c(0.8, 0.25))
})

ggs
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

Acoustic space projection

scale_acous_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median",
    "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
    "mindom", "maxdom", "dfrange", "modindx", "meanpeakf")])

urfmod <- randomForest(x = scale_acous_param, importance = TRUE, proximity = TRUE,
    ntree = 10000)

saveRDS(urfmod, "./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

mds_urf_prox <- cmdscale(1 - urfmod$proximity, k = 2)

saveRDS(mds_urf_prox, "./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

Random forest

urfmod <- readRDS("./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

mds_urf_prox <- readRDS("./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

sels$MDS1 <- mds_urf_prox[, 1]
sels$MDS2 <- mds_urf_prox[, 2]

print(1 - sum(urfmod$votes[, 1] > urfmod$votes[, 2])/nrow(urfmod$votes))
## [1] 0.005058366

Proportion of real data classified as synthetic data by the unsupervised random forest: 0.005058

ggplot(sels, aes(x = MDS1, y = MDS2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
    scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 30) + theme(legend.position = c(0.62,
    0.3)) + guides(colour = guide_legend(override.aes = list(size = 10)))

PCA

pca <- prcomp(x = sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
    "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
    "maxdom", "dfrange", "modindx", "meanpeakf")], scale. = TRUE)

Y <- as.data.frame(pca$x[, 1:2])
names(Y) <- c("PC1", "PC2")

sels <- data.frame(sels, Y)

ggplot(sels, aes(x = PC1, y = PC2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
    scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) + theme(legend.position = c(0.9,
    0.8)) + guides(colour = guide_legend(override.aes = list(size = 10)))

t-SNE

scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
    "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
    "maxdom", "dfrange", "modindx", "meanpeakf")])

tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE, max_iter = 5000)

saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")

Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")

sels <- data.frame(sels, Y)

ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(round))) + geom_point() +
    labs(color = "Round") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
    theme(legend.position = c(0.955, 0.25)) + guides(colour = guide_legend(override.aes = list(size = 10)))

sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress", "High Stress"))

ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) + geom_point() +
    labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
    theme(legend.position = c(0.88, 0.2)) + guides(colour = guide_legend(override.aes = list(size = 10)))

Acoustic space by individual (each point) for the 3 dimension reduction methods

  • Only individuals with at least 20 calls
  • Using rarefaction with n = the minimum sample size across all individuals
tab <- table(sels$ID)

Y <- sels[sels$ID %in% names(tab[tab > 20]), ]

pca_areas <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID",
    parallel = 4, pb = TRUE, type = "density")

rf_mds_areas <- rarefact_space_size(X = Y, dimensions = c("MDS1", "MDS2"), group = "ID",
    parallel = 4, pb = TRUE, type = "density")

tsne_areas <- rarefact_space_size(X = Y, dimensions = c("TSNE1", "TSNE2"), group = "ID",
    parallel = 4, pb = TRUE, type = "density")

areas <- data.frame(pca_areas[, -ncol(pca_areas)], pca.area = pca_areas$mean.size,
    rf.mds.area = rf_mds_areas$mean.size, tsne.area = tsne_areas$mean.size)

saveRDS(areas, "./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")
areas <- readRDS("./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")

ggpairs(areas, columns = which(names(areas) %in% c("pca.area", "rf.mds.area", "tsne.area")),
    columnLabels = c("PCA", "Rand.for. MDS", "t-SNE")) + theme_minimal(base_size = 30)

Individual acoustic space across weeks

  • PCA and t-SNE used for dimensionality reduction
  • Only weeks with at least 6 calls were included

PCA

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")

tab <- table(sels$ID.week)

# sort(tab)

Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]

# run with rarefaction
areas_by_week_raref <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
    parallel = 1, pb = TRUE, type = "density")

# run without rarefaction
areas_by_week <- space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
    parallel = 1, pb = TRUE, type = "density", extent = 2)

areas_by_week$raref.area <- areas_by_week_raref$mean.size

# add metadata
areas_by_week$ID <- sapply(strsplit(areas_by_week$group, "-"), "[[", 1)

areas_by_week$week <- factor(sapply(strsplit(areas_by_week$group, "-"), "[[", 2),
    levels = 1:5)

areas_by_week$treatment <- sapply(1:nrow(areas_by_week), function(x) metadat$Treatment[metadat$Bird.ID ==
    areas_by_week$ID[x]][1])

areas_by_week$round <- sapply(1:nrow(areas_by_week), function(x) metadat$Round[metadat$Bird.ID ==
    areas_by_week$ID[x]][1])

areas_by_week$round <- paste("Round", areas_by_week$round)

call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

areas_by_week$treatment <- factor(areas_by_week$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

saveRDS(areas_by_week, "./data/processed/acoustic_space_area_by_individual_and_week.RDS")

With rarefaction

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")

ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

With rarefaction and compared to week 1

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")

areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {

    X <- areas_by_week[areas_by_week$ID == x, ]
    X$raref.area <- X$raref.area - X$raref.area[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))

ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction

ggplot(data = areas_by_week[areas_by_week$week != 5, ], aes(x = week, y = size, group = ID,
    col = treatment)) + geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction and compared to week 1

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")

areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {

    X <- areas_by_week[areas_by_week$ID == x, ]
    X$size <- X$size - X$size[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))

ggplot(data = areas_by_week, aes(x = week, y = size, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Acoustic space area by number of calls per individual and week

ggplot(data = areas_by_week, aes(x = n, y = size, col = treatment)) + geom_point() +
    # geom_smooth(size = 1.2, method = 'lm') +
scale_color_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~round, ncol = 2, scales = "free_x") +
    labs(y = "Acoustic space area", x = "Number of calls (n)") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.25))

  • Acoustic space area by number of calls per individual and week by treatment
  • Must have at least 2 calls
ggplot(data = areas_by_week, aes(x = n, y = log(size + 10), col = treatment)) + geom_point() +
    geom_smooth(size = 1.2, method = "lm") + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    # facet_wrap(~ round, ncol = 2, scales = 'free_x') +
labs(y = "log(acoustic space area)", x = "Number of calls (n)") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.7))

Save acoustic space plots

trait_space_plot(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
    "TSNE2"), point.alpha = 0.2, pch = 1)
#
trait_space_plot(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
    "TSNE2"), background.indices = which(sels$ID != "395WBHM" & sels$record.group ==
    sels$record.group[sels$ID == "395WBHM"][1]), point.alpha = 0.2, pch = 1)
# Y <- sels[sels$ID == '395WBHM',] Y <- sels[sels$ID %in% c('395WBHM',
# '394WBHM', '360YGHM', '385YBHM'), ]

# trait_overlap(Y, dimensions = c('TSNE1', 'TSNE2'), group = 'ID', type =
# 'density', overlap.type = 'assymetric')


# X = Y dimensions = c('TSNE1', 'TSNE2') level <- '395WBHM' basecex <- 1 group
# <- 'ID' indices <- which(X[, group] == level)

trait_space_plot(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
    "TSNE2"), background.indices = which(sels$ID != "395WBHM" & sels$record.group ==
    sels$record.group[sels$ID == "395WBHM"][1]), point.alpha = 0.2, pch = 1)

# with points
out <- pblapply(unique(sels$ID), function(x) {

    tiff(filename = file.path("./images", paste(x, sels$treatment[sels$ID == x][1],
        "v2.tiff", sep = "-")), res = 120, width = 800, height = 800)

    par(mfrow = c(2, 2))

    # W <- sels[sels$ID == x, ]

    for (i in 1:4) try(trait_space_plot(X = sels, dimensions = c("TSNE1", "TSNE2"),
        indices = which(sels$ID == x & sels$week == i), background.indices = which(sels$ID !=
            x & sels$week == i & sels$record.group == sels$record.group[sels$ID ==
            x][1]), title = paste(x, "week", i), basecex = 1, labels = c(x, "group"),
        legend = FALSE), silent = TRUE)

    dev.off()
})

TSNE

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")

tab <- table(sels$ID.week)

# sort(tab)

Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]

# run with rarefaction
areas_by_week_raref <- rarefact_space_size(X = Y, dimensions = c("TSNE1", "TSNE2"),
    group = "ID.week", parallel = 1, pb = TRUE, type = "density")

# run without rarefaction
areas_by_week <- space_size(X = Y, dimensions = c("TSNE1", "TSNE2"), group = "ID.week",
    parallel = 4, pb = TRUE, type = "density")

areas_by_week$raref.area <- areas_by_week_raref$mean.size

# add metadata
areas_by_week$ID <- sapply(strsplit(areas_by_week$group, "-"), "[[", 1)

areas_by_week$week <- factor(sapply(strsplit(areas_by_week$group, "-"), "[[", 2),
    levels = 1:5)

areas_by_week$treatment <- sapply(1:nrow(areas_by_week), function(x) metadat$Treatment[metadat$Bird.ID ==
    areas_by_week$ID[x]][1])

areas_by_week$round <- sapply(1:nrow(areas_by_week), function(x) metadat$Round[metadat$Bird.ID ==
    areas_by_week$ID[x]][1])

areas_by_week$round <- paste("Round", areas_by_week$round)

call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

areas_by_week$treatment <- factor(areas_by_week$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

saveRDS(areas_by_week, "./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")

With rarefaction

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")

ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

With rarefaction and compared to week 1

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")

areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {

    X <- areas_by_week[areas_by_week$ID == x, ]
    X$raref.area <- X$raref.area - X$raref.area[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))

ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction

ggplot(data = areas_by_week[areas_by_week$week != 5, ], aes(x = week, y = size, group = ID,
    col = treatment)) + geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction and compared to week 1

areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {

    X <- areas_by_week[areas_by_week$ID == x, ]
    X$size <- X$size - X$size[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))

ggplot(data = areas_by_week, aes(x = week, y = size, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Acoustic space area by number of calls per individual and week

ggplot(data = areas_by_week, aes(x = n, y = size, col = treatment)) + geom_point() +
    # geom_smooth(size = 1.2, method = 'lm') +
scale_color_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~round, ncol = 2, scales = "free_x") +
    labs(y = "Acoustic space area", x = "Number of calls (n)") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.25))

  • Acoustic space area by number of calls per individual and week by treatment
  • Must have at least 2 calls
ggplot(data = areas_by_week, aes(x = n, y = log(size + 10), col = treatment)) + geom_point() +
    geom_smooth(size = 1.2, method = "lm") + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    # facet_wrap(~ round, ncol = 2, scales = 'free_x') +
labs(y = "log(acoustic space area)", x = "Number of calls (n)") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.7))

Acoustic space overlap

  • 2 methods: pairwise distance between observations and kernel density overlap

Individual overlap to its own group

PCA

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]

group_ovlp_l <- lapply(unique(Y$ID.week), function(x) {

    # group data
    X <- Y[Y$record.group == Y$record.group[Y$ID.week == x][1] & Y$week == Y$week[Y$ID.week ==
        x][1], ]

    # label all other individuals as group
    X$ID2 <- ifelse(X$ID.week == x, "focal", "group")

    # run with density
    ovlp <- if (min(table(X$ID2)) > 4)
        rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2",
            parallel = 1, pb = FALSE, outliers = 0.95, type = "mean.density.overlap") else matrix(rep(NA, 3), nrow = 1)

    dsts <- if (min(table(X$ID2)) > 1)
        rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2",
            parallel = 1, pb = FALSE, outliers = 0.95, type = "distance") else matrix(rep(NA, 3), nrow = 1)

    out <- data.frame(ID = Y$ID[Y$ID.week == x][1], group = Y$record.group[Y$ID.week ==
        x][1], week = Y$week[Y$ID.week == x][1], overlap.to.group = ovlp[1, 3], distance.to.group = dsts[1,
        3], treatment = Y$treatment[Y$ID.week == x][1], round = Y$round[Y$ID.week ==
        x][1])

    return(out)

})

group_ovlp <- do.call(rbind, group_ovlp_l)

saveRDS(group_ovlp, "./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")


group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID), function(x) {

    X <- group_ovlp[group_ovlp$ID == x, ]
    X$overlap.to.group <- X$overlap.to.group - X$overlap.to.group[X$week == min(as.numeric(as.character(X$week)))]
    X$distance.to.group <- X$distance.to.group - X$distance.to.group[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))


group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to group acoustic space",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to group acoustic space",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")

tab <- table(sels$ID.week)

Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]

# run with rarefaction
overlap_by_week <- rarefact_space_similarity(X = Y, dimensions = c("PC1", "PC2"),
    group = "ID.week", parallel = 1, pb = TRUE, outliers = 0.95, type = "mean.density.overlap")

saveRDS(overlap_by_week, "./data/processed/acoustic_space_overlap_by_individual_and_week.RDS")

TSNE

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 6], ]

group_ovlp_l <- lapply(unique(Y$ID.week), function(x) {

    # group data
    X <- Y[Y$record.group == Y$record.group[Y$ID.week == x][1] & Y$week == Y$week[Y$ID.week ==
        x][1], ]

    # label all other individuals as group
    X$ID2 <- ifelse(X$ID.week == x, "focal", "group")

    # run with density
    ovlp <- if (min(table(X$ID2)) >= 6)
        rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"), group = "ID2",
            parallel = 1, pb = FALSE, type = "mean.density.overlap")[, 1:3] else matrix(rep(NA, 3), nrow = 1)

    dsts <- if (min(table(X$ID2)) >= 2)
        rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"), group = "ID2",
            parallel = 1, pb = FALSE, type = "distance")[, 1:3] else matrix(rep(NA, 3), nrow = 1)

    out <- data.frame(ID = X$ID[X$ID.week == x][1], group = X$record.group[X$ID.week ==
        x][1], week = X$week[X$ID.week == x][1], overlap.to.group = ovlp[1, 3], distance.to.group = dsts[1,
        3], treatment = X$treatment[X$ID.week == x][1], round = X$round[X$ID.week ==
        x][1])

    return(out)

})

group_ovlp <- do.call(rbind, group_ovlp_l)


saveRDS(group_ovlp, "./data/processed/acoustic_space_density_overlap_to_group_by_week_TSNE.RDS")
group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week_TSNE.RDS")


# comparing to week 1 group_ovlp <- do.call(rbind,
# lapply(unique(group_ovlp$ID), function(x){ X <- group_ovlp[group_ovlp$ID ==
# x, ] X$overlap.to.group <- X$overlap.to.group / X$overlap.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] X$distance.to.group <-
# X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))


group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to group acoustic space",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

Change in overlap with self over time

PCA

tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 6], ]

indiv_ovlp_l <- lapply(unique(Y$ID), function(x) {

    # group data
    X <- Y[Y$ID == x, ]
    tab_week <- table(X$week)
    X <- X[X$week %in% names(tab_week)[tab_week >= 6], ]

    if (any(X$week == 1)) {
        # run with density
        ovlp <- try(rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"),
            group = "week", parallel = 1, pb = FALSE, type = "mean.density.overlap"),
            silent = TRUE)

        dsts <- try(rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"),
            group = "week", parallel = 1, pb = FALSE, type = "distance"), silent = TRUE)

        out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = if (is(ovlp,
            "try-error"))
            NA else ovlp$group.2[ovlp$group.1 == 1], overlap.to.first.week = if (is(ovlp,
            "try-error"))
            NA else ovlp$mean.overlap[ovlp$group.1 == 1], distance.to.first.week = if (is(dsts,
            "try-error"))
            NA else dsts$mean.distance[dsts$group.1 == 1], treatment = X$treatment[X$ID ==
            x][1], round = X$round[X$ID == x][1])
    } else out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = NA,
        overlap.to.first.week = NA, distance.to.first.week = NA, treatment = X$treatment[X$ID ==
            x][1], round = X$round[X$ID == x][1])
    return(out)

})

indiv_ovlp <- do.call(rbind, indiv_ovlp_l)

saveRDS(indiv_ovlp, "./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")
indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")

indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))
indiv_ovlp <- indiv_ovlp[!is.na(indiv_ovlp$group), ]

ggplot(data = indiv_ovlp, aes(x = week, y = overlap.to.first.week, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to acoustic space from first week",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

ggplot(data = indiv_ovlp, aes(x = week, y = distance.to.first.week, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to acoustic space from first week",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

agg_dist <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_dist$se <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
    3]
agg_dist$Week <- 1:4

# ggplot(data = agg_dist, aes(x = Week, y = distance.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=distance.to.first.week - se, ymax =
# distance.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Distance to self first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')


agg_ovlp <- aggregate(overlap.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_ovlp$se <- aggregate(overlap.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
    3]
agg_ovlp$Week <- 1:4

# ggplot(data = agg_ovlp, aes(x = Week, y = overlap.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=overlap.to.first.week - se, ymax =
# overlap.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Overlap to \nself first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')

TSNE

tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]

indiv_ovlp_l <- lapply(unique(Y$ID), function(x) {

    # group data
    X <- Y[Y$ID == x, ]
    tab_week <- table(X$week)
    X <- X[X$week %in% names(tab_week)[tab_week >= 6], ]

    if (any(X$week == 1)) {
        # run with density
        ovlp <- try(rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"),
            group = "week", parallel = 1, pb = FALSE, type = "mean.density.overlap"),
            silent = TRUE)

        dsts <- try(rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"),
            group = "week", parallel = 1, pb = FALSE, type = "distance"), silent = TRUE)

        out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = if (is(ovlp,
            "try-error"))
            NA else ovlp$group.2[ovlp$group.1 == 1], overlap.to.first.week = if (is(ovlp,
            "try-error"))
            NA else ovlp$mean.overlap[ovlp$group.1 == 1], distance.to.first.week = if (is(dsts,
            "try-error"))
            NA else dsts$mean.distance[dsts$group.1 == 1], treatment = X$treatment[X$ID ==
            x][1], round = X$round[X$ID == x][1])
    } else out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = NA,
        overlap.to.first.week = NA, distance.to.first.week = NA, treatment = X$treatment[X$ID ==
            x][1], round = X$round[X$ID == x][1])
    return(out)

})

indiv_ovlp <- do.call(rbind, indiv_ovlp_l)

saveRDS(indiv_ovlp, "./data/processed/acoustic_space_density_overlap_to_first_week_by_individual_TSNE.RDS")
indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual_TSNE.RDS")

indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))
indiv_ovlp <- indiv_ovlp[!is.na(indiv_ovlp$group), ]

ggplot(data = indiv_ovlp, aes(x = week, y = overlap.to.first.week, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to acoustic space from first week",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

ggplot(data = indiv_ovlp, aes(x = week, y = distance.to.first.week, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to acoustic space from first week",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

agg_dist <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_dist$se <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
    3]
# agg_stress$Week <- 1:4


# ggplot(data = agg_dist, aes(x = week, y = distance.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=distance.to.first.week - se, ymax =
# distance.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Distance to first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')

Statistical analyses

  • on t-SNE dimensions
  • predictors were mean-centered
  • link functions: binomial distribution for overlap and gaussian for area and call count
  • cort.stress and cort.baseline were dropped as they are highly correlated to stress response
  • Not a lot of data makes having many predictors complicated
  • No random factor (only 1 row per individual)
  • No p values printed. Look at whether credibility intervals (l-95% CI & u-95% CI) include 0 to assess statistical significance

Overlap to self (week 1) as response

overlap.to.first.week ~ treatment + stress.response

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
indiv_ovlp$ID.week <- paste(indiv_ovlp$ID, indiv_ovlp$week, sep = "-")

# addd stress parameters
indiv_ovlp$cort.baseline <- sapply(indiv_ovlp$ID.week, function(x) sels$cort.baseline[sels$ID.week ==
    x][1])
indiv_ovlp$cort.stress <- sapply(indiv_ovlp$ID.week, function(x) sels$cort.stress[sels$ID.week ==
    x][1])
indiv_ovlp$breath.count <- sapply(indiv_ovlp$ID.week, function(x) sels$breath.count[sels$ID.week ==
    x][1])
indiv_ovlp$stress.response <- indiv_ovlp$cort.stress - indiv_ovlp$cort.baseline

# mean centering
indiv_ovlp$cort.baseline <- mean(indiv_ovlp$cort.baseline, na.rm = TRUE) - indiv_ovlp$cort.baseline
indiv_ovlp$cort.stress <- mean(indiv_ovlp$cort.stress, na.rm = TRUE) - indiv_ovlp$cort.stress
indiv_ovlp$breath.count <- mean(indiv_ovlp$breath.count, na.rm = TRUE) - indiv_ovlp$breath.count
indiv_ovlp$stress.response <- mean(indiv_ovlp$stress.response, na.rm = TRUE) - indiv_ovlp$stress.response


# only individuals present in week 1 indiv_ovlp <- indiv_ovlp[indiv_ovlp$ID
# %in% indiv_ovlp$ID[indiv_ovlp$week == 1], ]

# indiv_ovlp <- indiv_ovlp[complete.cases(indiv_ovlp), ]

indiv_ovlp$overlap.to.first.week[indiv_ovlp$overlap.to.first.week == 0] <- 1e-04
indiv_ovlp$overlap.to.first.week[indiv_ovlp$overlap.to.first.week == 1] <- 0.9999

week 2:

## week 2

overlap_week_2 <- brm(iter = 5000, silent = 2, overlap.to.first.week ~ treatment +
    stress.response, data = indiv_ovlp[indiv_ovlp$week == 2, ], family = Beta(),
    control = list(adapt_delta = 0.9), chains = 1)

saveRDS(overlap_week_2, "./data/processed/brms_model_overlap_week_2.RDS")
overlap_week_2 <- readRDS("./data/processed/brms_model_overlap_week_2.RDS")

summary_brm_model(overlap_week_2)
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
treatmentMediumStress 0.065 1 1527.558 5000 1 -0.723 0.798
treatmentHighStress 0.618 1 1453.076 5000 1 -0.131 1.334
stress.response -0.006 1 2615.287 5000 1 -0.026 0.013

week 4:

## week 4
overlap_week_4 <- brm(iter = 5000, silent = 2, overlap.to.first.week ~ treatment +
    stress.response, data = indiv_ovlp[indiv_ovlp$week == 4, ], family = Beta(),
    control = list(adapt_delta = 0.9, max_treedepth = 12), chains = 1)

saveRDS(overlap_week_4, "./data/processed/brms_model_overlap_week_4.RDS")
overlap_week_4 <- readRDS("./data/processed/brms_model_overlap_week_4.RDS")

summary_brm_model(overlap_week_4)
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
treatmentMediumStress -0.702 1 2257.376 5000 1 -1.971 0.613
treatmentHighStress 0.454 1 2095.748 5000 1 -0.639 1.556
stress.response 0.042 1 2357.324 5000 1 0.001 0.082

Acoustic space change (compare to week 1)

area ~ treatment + stress.response

areas_by_week$ID.week <- paste(areas_by_week$ID, areas_by_week$week, sep = "-")


# addd stress parameters
areas_by_week$cort.baseline <- sapply(areas_by_week$ID.week, function(x) sels$cort.baseline[sels$ID.week ==
    x][1])
areas_by_week$cort.stress <- sapply(areas_by_week$ID.week, function(x) sels$cort.stress[sels$ID.week ==
    x][1])
areas_by_week$breath.count <- sapply(areas_by_week$ID.week, function(x) sels$breath.count[sels$ID.week ==
    x][1])
areas_by_week$stress.response <- areas_by_week$cort.stress - areas_by_week$cort.baseline

# mean centering
areas_by_week$cort.baseline <- mean(areas_by_week$cort.baseline, na.rm = TRUE) -
    areas_by_week$cort.baseline
areas_by_week$cort.stress <- mean(areas_by_week$cort.stress, na.rm = TRUE) - areas_by_week$cort.stress
areas_by_week$breath.count <- mean(areas_by_week$breath.count, na.rm = TRUE) - areas_by_week$breath.count
areas_by_week$stress.response <- mean(areas_by_week$stress.response, na.rm = TRUE) -
    areas_by_week$stress.response


# only individuals present in week 1 areas_by_week <-
# areas_by_week[areas_by_week$ID %in% areas_by_week$ID[areas_by_week$week ==
# 1], ]


# areas_by_week <- areas_by_week[complete.cases(areas_by_week), ]

week 2:

## week 2
area_week_2 <- brm(iter = 5000, silent = 2, size ~ treatment + stress.response, data = areas_by_week[areas_by_week$week ==
    2, ], control = list(adapt_delta = 0.9), chains = 1)

saveRDS(area_week_2, "./data/processed/brms_model_area_week_2.RDS")
area_week_2 <- readRDS("./data/processed/brms_model_area_week_2.RDS")

summary_brm_model(area_week_2)
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
treatmentMediumStress 0.175 1 1787.606 5000 1 -0.314 0.672
treatmentHighStress 0.126 1 1701.261 5000 1 -0.384 0.614
stress.response 0.012 1 2390.428 5000 1 -0.002 0.026

week 4:

## week 4
area_week_4 <- brm(iter = 5000, silent = 2, size ~ treatment + stress.response, data = areas_by_week[areas_by_week$week ==
    4, ], control = list(adapt_delta = 0.9), chains = 1)

saveRDS(area_week_4, "./data/processed/brms_model_area_week_4.RDS")
area_week_4 <- readRDS("./data/processed/brms_model_area_week_4.RDS")

summary_brm_model(area_week_4)
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
treatmentMediumStress 0.461 1.004 1131.119 5000 1 -0.335 1.255
treatmentHighStress 0.669 1.006 1362.654 5000 1 -0.054 1.394
stress.response 0.017 1.001 1806.867 5000 1 -0.013 0.045

Call count

week 2:

## week 2
count_week_2 <- brm(iter = 5000, silent = 2, selec ~ treatment + stress.response,
    data = call_count[call_count$week == 2, ], control = list(adapt_delta = 0.9),
    chains = 1)

saveRDS(count_week_2, "./data/processed/brms_model_count_week_2.RDS")
count_week_2 <- readRDS("./data/processed/brms_model_count_week_2.RDS")

summary_brm_model(count_week_2)
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
treatmentMediumStress 355.784 1.003 1074.230 5000 1 -127.180 804.368
treatmentHighStress 613.795 1.001 1050.304 5000 1 131.824 1086.057
stress.response 0.971 1 1396.041 5000 1 -9.032 10.761

week 4:

## week 4
count_week_4 <- brm(iter = 5000, silent = 2, selec ~ treatment + stress.response,
    data = call_count[call_count$week == 4, ], control = list(adapt_delta = 0.9),
    chains = 1)

saveRDS(count_week_4, "./data/processed/brms_model_count_week_4.RDS")
count_week_4 <- readRDS("./data/processed/brms_model_count_week_4.RDS")

summary_brm_model(count_week_4)
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
treatmentMediumStress -103.093 1 1001.866 5000 1 -335.787 123.953
treatmentHighStress 3.989 1.001 840.242 5000 1 -352.231 359.975
stress.response 2.278 1.001 897.594 5000 1 -1.867 5.900

 

Takeaways

  • Stress has an effect on number of calls: high-stress individuals produce more calls than low-stress individuals in week 2
  • Stress response has a positive effect on the overlap to week 1 (overlap to itself)

 


Proposed Models for Stress and Learning Pilot Data based on 9/28/21 discussion

Model set 1: effects of treatment

Predicted variables: - Physiology: Baseline CORT, stress response CORT (stress), difference baseline-stress response CORT (stress response), weight, breath rate - Learning: number of calls, total acoustic area, distance moved from week 1, overlap to self week 1, overlap self to group

General form: - Predicted ~ treatment + week + ind(rand) - Predicted ~ treatment + week + treat * week + ind(rand) - Predicted ~ round + week [to check for effect of round]

Note: Use week 1 as baseline, comparing weeks 2, 3, 4

agg_dat <- aggregate(selec ~ ID + week, data = sels, length)

# compare to week 1 agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) {
# baseline <- agg_dat$selec[agg_dat$week == 1 & agg_dat$ID == agg_dat$ID[x]] if
# (length(baseline) > 0) change <- agg_dat$selec[x] - baseline else change <-
# agg_dat$selec[x] return(change) } )

# without comparing to week 1
agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) agg_dat$selec[x])

agg_dat$selec <- NULL

agg_dat$baseline.CORT <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$cort.baseline[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$cort.baseline[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$stress.response <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$stress.response[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$stress.response[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})


agg_dat$stress.CORT <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$cort.stress[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$cort.stress[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$weight <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$weight[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$weight[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$breath.rate <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$breath.count[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$breath.count[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$acoustic.area <- sapply(1:nrow(agg_dat), function(x) {

    area <- areas_by_week$raref.area[areas_by_week$ID == agg_dat$ID[x] & areas_by_week$week ==
        agg_dat$week[x]]

    if (length(area) < 1)
        area <- NA

    return(area)
})


agg_dat$acoustic.distance <- sapply(1:nrow(agg_dat), function(x) {

    distance <- indiv_ovlp$distance.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
        indiv_ovlp$week == agg_dat$week[x]]

    if (length(distance) < 1)
        distance <- NA

    return(distance)
})

agg_dat$acoustic.overlap <- sapply(1:nrow(agg_dat), function(x) {

    overlap <- indiv_ovlp$overlap.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
        indiv_ovlp$week == agg_dat$week[x]]

    if (length(overlap) < 1)
        overlap <- NA

    return(overlap)
})

agg_dat$acoustic.overlap.to.group <- sapply(1:nrow(agg_dat), function(x) {

    overlap <- group_ovlp$overlap.to.group[group_ovlp$ID == agg_dat$ID[x] & group_ovlp$week ==
        agg_dat$week[x]]

    if (length(overlap) < 1)
        overlap <- NA

    return(overlap)
})


agg_dat$treatment <- sapply(1:nrow(agg_dat), function(x) unique(sels$treatment[sels$ID ==
    agg_dat$ID[x]]))

agg_dat$round <- sapply(1:nrow(agg_dat), function(x) unique(sels$round[sels$ID ==
    agg_dat$ID[x]]))

# remove week 1
agg_dat <- agg_dat[agg_dat$week != 1, ]


responses <- setdiff(names(agg_dat), c("ID", "week", "treatment", "round"))

predictors <- c("~ treatment + week + (1|ID)", "~ treatment * week + (1|ID)", "~ round + week + (1|ID)")

formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
vars_to_scale <- names(agg_dat)[sapply(agg_dat, is.numeric)]
vars_to_scale <- vars_to_scale[!vars_to_scale %in% c("week", "round")]

for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[, vars_to_scale])

treatment_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
        ]
    sub_dat

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 2), silent = TRUE)
    return(mod)
})


names(treatment_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(treatment_models, "./data/processed/treatment_models.RDS")

agg_dat$week <- as.factor(agg_dat$week)

treatment_models_week_factor <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
        ]

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 2), silent = TRUE)
    return(mod)
})

names(treatment_models_week_factor) <- paste(formulas$responses, formulas$predictors)

saveRDS(treatment_models_week_factor, "./data/processed/treatment_models_week_factor.RDS")

week as continuous

treatment_models <- readRDS("./data/processed/treatment_models.RDS")



for (x in treatment_models) {
    summary_brm_model(x)
}
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT -0.342 1.001 2032.331 5000 2 -0.623 -0.047
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT 0.057 1 3386.12 5000 2 -0.198 0.305
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT -0.08 1.001 2172.071 5000 2 -0.368 0.209
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT 0.018 1.001 3184.54 5000 2 -0.28 0.314
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT 0.062 1 2669.238 5000 2 -0.227 0.366
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response 0.267 1.001 1967.507 5000 2 -0.04 0.569
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response -0.14 1 3307.176 5000 2 -0.397 0.111
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response 0.067 1 2819.069 5000 2 -0.213 0.339
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response -0.089 1 2859.569 5000 2 -0.384 0.212
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response 0.001 1 3164.738 5000 2 -0.309 0.302

 

Takeaways

  • No significant nor marginal effects

 

Model set 2: individual stress and effects on learning

General form: - Learning variables ~ baseline CORT [Week 2 immediate stress and learning, week 3 chronic stress and learn]

agg_dat <- agg_dat[agg_dat$week != "4", ]

agg_dat$week <- as.factor(agg_dat$week)


responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
    "acoustic.overlap.to.group")

predictors <- c("~ baseline.CORT + week + (1|ID)", "~ stress.response + week + (1|ID)")

formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
treatment_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
        ]

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 1), silent = TRUE)
    return(mod)
})


names(treatment_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(treatment_models, "./data/processed/stress_on_learning_models.RDS")
treatment_models <- readRDS("./data/processed/stress_on_learning_models.RDS")



for (x in treatment_models) {
    summary_brm_model(x)
}
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT -0.192 1.002 997.111 5000 1 -0.546 0.177
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT -0.154 1 537.359 5000 1 -0.462 0.159
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT 0.007 1.015 256.643 5000 1 -0.248 0.274
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT 0.101 1 616.351 5000 1 -0.247 0.44
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response 0.299 1.001 1394.371 5000 1 -0.112 0.69
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response 0.013 1 556.089 5000 1 -0.25 0.287
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response -0.058 1.001 677.164 5000 1 -0.391 0.286
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response -0.138 1.002 869.593 5000 1 -0.49 0.209

 

Takeaways

  • No significant effects

 

Model set 3: individual stress and effects on learning on week 2 and 3 independently

Week 2 immediate stress and learning, week 3 chronic stress and learn

General form: - Learning variables ~ baseline CORT - Learning variables ~ difference baseline-stress CORT

Week 2

agg_dat_w2 <- agg_dat[agg_dat$week == "2", ]

responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
    "acoustic.overlap.to.group")

predictors <- c("~ baseline.CORT + (1|ID)", "~ stress.response + (1|ID)")

formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)

treatment_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat_w2[!is.na(agg_dat_w2[names(agg_dat_w2) == formulas$responses[x]]),
        ]

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 1), silent = TRUE)
    return(mod)
})


names(treatment_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(treatment_models, "./data/processed/stress_on_learning_models_week_2.RDS")
treatment_models <- readRDS("./data/processed/stress_on_learning_models_week_2.RDS")

for (x in treatment_models) {
    summary_brm_model(x)
}
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT -0.177 1 1201.316 5000 1 -0.501 0.164
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT 0.037 1.002 701.318 5000 1 -0.328 0.408
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT -0.138 1.012 135.412 5000 1 -0.426 0.133
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT -0.03 1.045 78.7 5000 1 -0.329 0.298
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT 0.11 1.004 1054.756 5000 1 -0.257 0.458
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response 0.261 1.001 831.794 5000 1 -0.072 0.586
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response -0.178 1 227.403 5000 1 -0.545 0.202
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response 0.013 1 1041.844 5000 1 -0.232 0.25
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response -0.045 1.002 546.832 5000 1 -0.381 0.267
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response -0.138 1.004 1319.857 5000 1 -0.529 0.261

Week 3

agg_dat_w3 <- agg_dat[agg_dat$week == "3", ]

responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
    "acoustic.overlap.to.group")

predictors <- c("~ baseline.CORT + (1|ID)", "~ stress.response + (1|ID)")

formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)

treatment_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat_w3[!is.na(agg_dat_w3[names(agg_dat_w3) == formulas$responses[x]]),
        ]

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 1), silent = TRUE)
    return(mod)
})


names(treatment_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(treatment_models, "./data/processed/stress_on_learning_models_week_3.RDS")
treatment_models <- readRDS("./data/processed/stress_on_learning_models_week_3.RDS")

for (x in treatment_models) {
    summary_brm_model(x)
}
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT -0.192 1.002 997.111 5000 1 -0.546 0.177
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT -0.154 1 537.359 5000 1 -0.462 0.159
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT 0.007 1.015 256.643 5000 1 -0.248 0.274
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
baseline.CORT 0.101 1 616.351 5000 1 -0.247 0.44
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response 0.299 1.001 1394.371 5000 1 -0.112 0.69
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response 0.013 1 556.089 5000 1 -0.25 0.287
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response -0.058 1.001 677.164 5000 1 -0.391 0.286
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
stress.response -0.138 1.002 869.593 5000 1 -0.49 0.209

 

Takeaways

  • No significant nor marginal effects

 


Session information

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] PhenotypeSpace_0.1.0  warbleR_1.1.26        NatureSounds_1.0.4   
##  [4] seewave_2.1.8         tuneR_1.3.3.1         ggridges_0.5.3       
##  [7] posterior_1.2.0       ggrepel_0.9.1         cowplot_1.1.1        
## [10] tidybayes_3.0.2       modelr_0.1.8          tidyr_1.2.0          
## [13] forcats_0.5.1         purrr_0.3.4           dplyr_1.0.8          
## [16] lme4_1.1-28           Matrix_1.3-4          brms_2.16.3          
## [19] Rcpp_1.0.8            vegan_2.5-7           lattice_0.20-44      
## [22] permute_0.9-7         raster_3.5-15         spatstat_2.3-3       
## [25] spatstat.linnet_2.3-2 spatstat.core_2.4-0   rpart_4.1-15         
## [28] nlme_3.1-152          spatstat.random_2.1-0 spatstat.geom_2.3-2  
## [31] spatstat.data_2.1-2   GGally_2.1.2          adehabitatHR_0.4.19  
## [34] adehabitatLT_0.3.25   CircStats_0.2-6       boot_1.3-28          
## [37] adehabitatMA_0.3.14   ade4_1.7-18           deldir_1.0-6         
## [40] sp_1.4-6              MASS_7.3-54           Rtsne_0.15           
## [43] randomForest_4.7-1    formatR_1.11          knitr_1.37           
## [46] kableExtra_1.3.4      ggplot2_3.3.5         viridis_0.6.2        
## [49] viridisLite_0.4.0     pbapply_1.5-0         readxl_1.3.1         
## [52] lubridate_1.8.0       remotes_2.4.2        
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2            tidyselect_1.1.2      fftw_1.0-6.1         
##   [4] htmlwidgets_1.5.4     grid_4.1.0            munsell_0.5.0        
##   [7] codetools_0.2-18      DT_0.20               miniUI_0.1.1.1       
##  [10] withr_2.4.3           Brobdingnag_1.2-7     colorspace_2.0-3     
##  [13] highr_0.9             rstudioapi_0.13       stats4_4.1.0         
##  [16] dtw_1.22-3            tensor_1.5            bayesplot_1.8.1      
##  [19] labeling_0.4.2        rstan_2.21.3          polyclip_1.10-0      
##  [22] farver_2.1.0          bridgesampling_1.1-2  coda_0.19-4          
##  [25] vctrs_0.3.8           generics_0.1.2        xfun_0.29            
##  [28] R6_2.5.1              markdown_1.1          HDInterval_0.2.2     
##  [31] gamm4_0.2-6           projpred_2.0.2        bitops_1.0-7         
##  [34] spatstat.utils_2.3-0  reshape_0.8.8         assertthat_0.2.1     
##  [37] promises_1.2.0.1      scales_1.1.1          gtable_0.3.0         
##  [40] processx_3.5.2        goftest_1.2-3         rlang_1.0.1          
##  [43] systemfonts_1.0.4     splines_4.1.0         broom_0.7.12         
##  [46] checkmate_2.0.0       inline_0.3.19         yaml_2.3.5           
##  [49] reshape2_1.4.4        abind_1.4-5           threejs_0.3.3        
##  [52] crosstalk_1.2.0       backports_1.4.1       httpuv_1.6.5         
##  [55] rsconnect_0.8.25      tensorA_0.36.2        tools_4.1.0          
##  [58] ellipsis_0.3.2        jquerylib_0.1.4       RColorBrewer_1.1-2   
##  [61] proxy_0.4-26          plyr_1.8.6            base64enc_0.1-3      
##  [64] RCurl_1.98-1.6        ps_1.6.0              prettyunits_1.1.1    
##  [67] zoo_1.8-9             cluster_2.1.2         magrittr_2.0.2       
##  [70] ggdist_3.1.0          colourpicker_1.1.1    mvtnorm_1.1-3        
##  [73] matrixStats_0.61.0    shinyjs_2.1.0         mime_0.12            
##  [76] evaluate_0.15         arrayhelpers_1.1-0    xtable_1.8-4         
##  [79] shinystan_2.5.0       gridExtra_2.3         rstantools_2.1.1     
##  [82] compiler_4.1.0        tibble_3.1.6          crayon_1.5.0         
##  [85] minqa_1.2.4           StanHeaders_2.21.0-7  htmltools_0.5.2      
##  [88] mgcv_1.8-36           later_1.3.0           RcppParallel_5.1.5   
##  [91] DBI_1.1.1             cli_3.2.0             parallel_4.1.0       
##  [94] igraph_1.2.11         pkgconfig_2.0.3       signal_0.7-7         
##  [97] terra_1.5-21          spatstat.sparse_2.1-0 xml2_1.3.3           
## [100] svUnit_1.0.6          dygraphs_1.1.1.6      svglite_2.1.0        
## [103] bslib_0.3.1           webshot_0.5.2         rvest_1.0.2          
## [106] stringr_1.4.0         distributional_0.3.0  callr_3.7.0          
## [109] digest_0.6.29         rmarkdown_2.11        cellranger_1.1.0     
## [112] shiny_1.7.1           gtools_3.9.2          rjson_0.2.21         
## [115] nloptr_2.0.0          lifecycle_1.0.1       jsonlite_1.8.0       
## [118] fansi_1.0.2           pillar_1.7.0          loo_2.4.1            
## [121] fastmap_1.1.0         httr_1.4.2            pkgbuild_1.3.1       
## [124] glue_1.6.1            xts_0.12.1            shinythemes_1.2.0    
## [127] stringi_1.7.6         sass_0.4.0